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I. INTRODUCTION 

In this paper we develop a general multi-mode Gaus- 
sian representation for a density matrix of bosons. As 
well as classical phase-space variables like (x, p) , the rep- 
resentation utilizes a dynamical space of quantum uncer- 
tainties or co-variances. The extended phase space ac- 
commodates more efficiently the content of a quantum 
state, and allows the physics of many kinds of problems 
to be incorporated into the basis itself. The Gaussian 
expansion technique unifies and greatly extends all the 
previous Gaussian-like phase-space representations used 
for bosons, including the Wigner, Q, P, positive-P, and 
squeezed-state expansions. The operator basis also in- 
cludes non-Hermitian Gaussian operators, which are not 
density matrices themselves, but can form part of a prob- 
abilistic expansion of a physical density matrix. Unlike 
previous approaches, any initial state is found to evolve 
with a deterministic time-evolution under any quadratic 
Hamiltonian or master equation. 

The complexity of many-body quantum physics is 
manifest in the enormity of the Hilbert space of systems 
with even modest numbers of particles. This complexity 
makes it prohibitively difficult to simulate quantum dy- 
namics with exact precision: no digital computer is large 
enough to store the dynamically evolving state. However, 
if a finite precision is permitted, then quantum dynam- 
ical calculations are possible, through what are known 
as phase-space methods. These methods represent the 
evolving quantum state as probability distributions on 
some suitable phase space, which can be sampled via 
stochastic techniques. The mapping to phase space can 
be made to be exact. Thus the precision of the final re- 
sult is limited only by sampling error, which can usually 
be reliably estimated and which can be reduced by an 
increased number of stochastic paths. 

Arbitrary quantum mechanical evolution cannot be 
represented (even probabilistically) on a phase space as 
is usually defined. Thus the extended phase spaces em- 
ployed here are a generalization of conventional phase 
spaces in several ways. First, it is a quantum phase space, 
in which points can correspond to states with intrinsic 
uncertainty. Heisenberg's uncertainty relations can thus 



be satisfied in this way, or more generally by consider- 
ing genuine probability distributions over phase space, 
to be sampled stochastically. Second, the phase space is 
of double dimension, where classically real variables, such 
as X and p, now range over the complex plane. This al- 
lows arbitrary quantum evolution to be sampled stochas- 
tically. Third, stochastic gauge variables are included. 
These arbitrary quantities do not affect the physical re- 
sults, but they can be used to overcome problems in the 
stochastic sampling. Fourth, the phase space includes the 
set of second-order moments, or covariances. A phase 
space that is enlarged in this way is able to accommo- 
date more information about a general quantum state in 
a single point. In particular, any state (pure or mixed) 
with Gaussian statistics can be represented as a single 
point in this phase space. 

The Gaussian representation provides a link between 
phase-space methods and approximate methods used in 
many-body theory, which frequently treat normal and 
anomalous correlations or Green's functions as dynami- 
cal objects As well as being applicable to quantum 
optics and quantum information, a strong motivation for 
this representation is the striking experimental observa- 
tion of BEG in ultra-cold atomic systems 14j. Already 
the term 'atom laser' is widely used, and experimen- 
tal observation of quantum statistics in these systems 
is underway. Yet there is a problem in using previous 
quantum optics formalisms to calculate coherence prop- 
erties in atom optics: interactions are generally much 
stronger with atoms than they are with photons, relative 
to the damping rate. The consequence of this is that 
one must anticipate larger departures from 'semiclassi- 
cal', coherent-state behaviour in atomic systems. 

The present paper includes these nonclassical and in- 
coherent effects at the level of the basis for the operator 
representation itself. The purpose of employing a Gaus- 
sian basis set is not only to enlarge the parameter set (to 
hold more information about the quantum state), but 
also to include in the basis states that are a close match 
to the actual states that are likely to occur in interesting 
systems, such as dilute gases. The pay-off for increas- 
ing the parameter set is more efficient sampling of the 
dynamically evolving or equilibrium states of many-body 
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systems. 

The idea of coherent states as a quasi-classical basis for 
quantum mechanics originated with SchrodingerQ. Sub- 
sequently, Wigner 2] introduced a distribution for quan- 
tum density matrices. This method was a phase-space 
mapping with classical dimensions and employed a sym- 
metrically ordered operator correspondence principle. 
Later developments included the antinormally ordered Q 
distribution 3] , a normally ordered expansion called the 
diagonal P distribution 4] , methods that interpolate be- 
tween these classical phase-space distributions^^, and 
diagonal squeezed-state representations 0. Each of these 
expansions either employs an explicit Gaussian density 
matrix basis or is related to one that does by convolu- 
tion. They are suitable for phase-space representations 
of quantum states because of the overcompleteness of the 
set of coherent states on which they are based. 

Arbitrary pure states of bosons with a Gaussian wave- 
function or Wigner representation are often called the 
squeezed states IS]. These are a superset of the coherent 
states, and were investigated by BogoliubovT^ to ap- 
proximately represent the ground state of an interacting 
Bose-Einstein condensate - as well as in much recent work 
in quantum optics T^. Diagonal expansions analogous to 
the diagonal P-rcpresentation have been introduced us- 
ing a basis of squeezed-state projectors, typically with a 
fixed squeezing parameter^- However, these have not 
generally resulted in useful dynamical applications, as 
they do not overcome the problems inherent in using a 
diagonal basis, as we discuss below. 

In operator representations, one must utilize a com- 
plete basis in the Hilbert space of density operators, 
rather than in the Hilbert space of pure states. Ther- 
mal density matrices, for example, are not pure states, 
but do have a Gaussian P representation and Wigner 
function. To include all three types of commonly used 
Gaussian states - the coherent, squeezed and thermal 
states - one can define a Gaussian state as a den- 
sity matrix having a Gaussian positive-P or Wigner 
representation [l^ . This definition also includes dis- 
placed and squeezed thermal states. Gaussian states 
have been investigated extensively in quantum informa- 
tion and quantum entanglement [19] . It has been shown 
that an initial Gaussian state will remain Gaussian under 
linear evolution[25j. 

However, the Gaussian density matrices that corre- 
spond to physical states do not by themselves form a 
complete basis for the time-evolution of all quantum den- 
sity matrices. This problem, inherent in all diagonal 
expansions, is related to known issues in constructing 
quantum-classical correspondences!^, and is caused by 
the non-positive-definite nature of the local propagator 
in a classical phase space. It is manifest in the fact that 
there is generally no equivalent Fokker-Planck equation 
(with a positive-definite diffusion matrix) that generates 
the quantum time-evolution, and hence no correspond- 
ing stochastic differential equation that can be efficiently 
simulated numerically. This difficulty occurs in nearly all 



cases except free fields, and represents a substantial limi- 
tation in the use of these diagonal expansion methods for 
exact simulation of the quantum dynamics of interacting 
systems. 

These problems can be solved by use of nonclassical 
phase spaces, which correspond to expansions in non- 
Hermitian bases of operators (rather than just physi- 
cal density matrices). One established example is the 
non-diagonal positive-PQ representation. The non- 
Hermitian basis in this case generates a representation 
with a positive propagator, which allows the use of 
stochastic methods to sample the quantum dynamics. 
By extending the expansion to include a stochastic gauge 
freedom in these evolution equations, one can select the 
most compact possible time-evolution equation^, 0, 
IT^ . With an appropriate gauge choice, this method is 
exact for a large class of nonlinear Hamiltonians, since it 
eliminates boundary terms that can otherwise arise [2 ij. 
The general Gaussian representation used here also in- 
cludes these features, and extends them to allow treat- 
ment of any Hamiltonian or master equation with up to 
fourth-order polynomial terms. 

Other methods of theoretical physics that have compa- 
rable goals are the path-integral techniques of quantum 
field theory [23, H^, and density functional methods 
that are widely used to treat atomic and molecular sys- 
tems. The first of these is exact in principle, but is al- 
most exclusively used in imaginary time calculations of 
canonical ensembles or ground states due to the notori- 
ous phase problem. The second method has similarities 
with our approach in that it also utilizes a density as 
we do. However, density functionals are normally com- 
bined with approximations like the local density approx- 
imation. Gaussian representation methods have the ad- 
vantage that they can treat both real and imaginary time 
evolution. In addition, the technique is exact in principle, 
provided boundary terms vanish on partial integration. 

In section m we define general Gaussian operators for 
a density operator expansion and introduce a compact 
notation for these operators, either in terms of mode op- 
erators or quantum fields. In section 11111 we calculate 
the moments of the general Gaussian representation, re- 
lating them to physical quantities as well as to the mo- 
ments of previous representations. Section llVI gives the 
necessary identities that enable first-principles quantum 
calculations with these representations. 

Equation 14.15|l summarizes the relevant operator 
mappings, and constitutes a key result of the paper. 

We give a number of examples in section [3 of specific 
pure and mixed states (and their non-Hermitian gener- 
alizations) that are included in the basis, and we give 
simplified versions of important identities for these cases. 
Section lVIl describes how the Gaussian representation can 
be used to deal with evolution in either real or imaginary 
time. In particular, we show how it can be used to solve 
exactly any master equation that is quadratic in annihi- 
lation and creation operators. Some useful normalization 
integrals and reordering identities for the Gaussian oper- 
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Figure 1: The density-operator expansion in Eq. H2.1^ can be 
interpreted as a convolution of the probability distribution P 
with the underlying distribution of the basis. The uncertainty 
or spread of the physical state, indicated by the variance ap, 
is shared between the distribution variance ap and the basis 
variance ua. 



ators are proved in the Appendices. 

In a subsequent paper, we will apply these methods to 
systems with nonlinear evolution. 



II. THE GAUSSIAN REPRESENTATION 

The representations that give exact mappings between 
operator equations and stochastic equations - an essen- 
tial step toward representing operator dynamics in large 
Hilbert spaces - are stochastic gauge expansions [ifl UM, 
IT^ on a nonclassical phase space. Here, the generic ex- 
pansion is written down in terms of a complete set of 
operators that are typically non-Hermitian. This leads 
to the typical form: 



P(A,i)A(A)dA 



(2.1) 



where P( A , t) is a probability distribution, A is a suitable 
basis for the class of density matrices being considered, 
and d A is the integration measure for the corresponding 
generalized phase-space coordinate A . See Fig. ^ for a 
conceptual illustration of this expansion. 

In phase-space methods, it is the distribution P that is 
sampled stochastically. Therefore if the basis resembles 
the typical physical states of a system, the sampling er- 
ror will be minimized, and if the state coincides exactly 
with an element of the basis, then the distribution will 
be a delta function, with consequently no sampling error. 
A Wigner or Q-function basis, for example, generates a 
broad distribution even for minimum uncertainty states. 
A general Gaussian basis, on the other hand, can gener- 
ate a delta-function distribution not only for any mini- 
mum uncertainty state, but also for the ground states of 
non-interacting finite-temperature systems. 



A. Gaussian operator basis 

In this paper, we define the operator basis A to be 
the most general Gaussian operator basis. The motiva- 
tion for using the most general possible basis set is that 
when the basis set members nearly match the states of in- 
terest, the resulting distributions are more compact and 
have lower sampling errors in a Monte-Carlo or stochas- 
tic calculation. In addition, a larger basis allows more 
choice of mappings, so that lower-order differential cor- 
respondences can be utilized. In some cases, a large basis 
set can increase computational memory requirements, as 
more parameters are needed. This disadvantage is out- 
weighed when there is a substantial decrease in sampling 
error, due to the use of a more physically appropriate 
basis. By choosing a general Gaussian operator basis, 
rather than just a basis of Gaussian density matrices, 
one has the additional advantage of a complete represen- 
tation for all non-Gaussian density matrices as well. 

If a is a column vector of M bosonic annihilation op- 
erators, and the corresponding row vector of creation 
operators, their commutation relations are: 



ak,a 



(2.2) 



We define a Gaussian operator as an exponential of an 
arbitrary quadratic form in annihilation and creation op- 
erators (or, equivalently, a quadratic form in position and 
momentum operators). 

The simplest way to achieve this is to introduce ex- 
tended 2M-vectors of c-numbers and operators: a = 
( a, (q:+)-^), and a = (2, (2^)^), with adjoints defined 
as a"*" = (Q;+,a-^) and = {a'',a'^), together with a 
relative operator displacement of 

/ 2i \ 
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(2.3) 



These extended vectors are indexed where necessary with 
Greek indices: fJ. — I, 2M. 

A general Gaussian operator is now an exponential of a 
general quadratic form in the 2M-vector mode operator 
So:. For algebraic reasons, it is useful to employ normal 
ordering, and to introduce a compact notation using a 
generalized covariance ct: 




Here the normalization factor involving y | ct | is intro- 
duced to simplify identities that occur later, and plays a 
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very similar role to the exactly analogous normalization 
factor that occurs in the classical Gaussian distribution 
of probability theory. The 2M x 2M covariance matrix 
is conveniently parameterized in terms of M x M subma- 
trices as: 



I 



n 



(2.5) 



where n is a complex M x M matrix and m, m+ are two 
independent symmetric complex M x M matrices. 

With this choice, the covariance has a type of gener- 
alized Hermitian symmetry in which ct^j^ = cr^^M,fj.+AU 
provided we interpret the matrix indices as cyclic in the 
sense that i/ ^ ly + 2M. This can also be written as 
g_ = CT+jWith the definition that: 



The complex vectors a and cc"*" give the generalized 
coherent amplitudes for each mode: a defines the ampli- 
tudes of annihilation operators a , while its 'conjugate' 
a"*" defines the amplitudes of the creation operators S^. 
The matrix n gives the number, or normal, correlations 
between each pair of modes. The squeezing, or anoma- 
lous, correlations between each pair of modes are given 
by m and m+: the matrix m defines the correlations of 
annihilation operators, while its 'conjugate' m+ defines 
the correlations of the creation operators. These physi- 
cal interpretations of the phase-space variables are sup- 
ported by the results of Section ITTll where we rigorously 
establish the connection of the phase-space variables to 
physical quantities. 

In general, apart from the complex amplitude fl, the 
total number of complex parameters needed to specify 
the normalized Af-mode Gaussian operator is: 



a b 
c d 



d c 
b a 



(2.6) 



This definition implies that we intend the superscript 
to define an operation on the covariance matrix which 
is equivalent to Hermitian conjugation of the underlying 
operators. If the Gaussian operator is in fact an Her- 
mitian operator, then so is the corresponding covariance 
matrix. In this case, the '"^' superscript is identical to 
ordinary Hermitian conjugation. The generalized Hermi- 
tian symmetry of the covariance means that all elements 
of the number correlation n appear twice, as do all ex- 
cept the diagonal elements of the squeezing correlations 
m, m+. 

The use of normal ordering allows simple operator 
identities to be obtained, but can easily be related to 
more commonly used unordered parameterizations. The 
Gaussian operators include as special cases the den- 
sity matrices of many useful and well-known physical 
states. For example, they include the thermal states of a 
Bose-Einstein distribution, the coherent states, and the 
squeezed states. They also include many more states 
than these, like the off-diagonal coherent state projectors 
used in the positive-P expansion, which are not density 
operators themselves, but can be used to expand density 
operators. The details are given in section Ivl 



B. Extended phase space 

The representation phase space is thus extended to 

X = (ri, a,a+,n, m, m+) . (2.7) 

The complex amplitude Q, which appears in the normal- 
isation, acts as a dynamical weight on different stochas- 
tic trajectories. It is useful in calculations in which the 
normalisation of the density matrix is not intrinsically 
preserved, such as canonical ensemble calculations, and 
also enables stochastic gauges to be included. 



p = M{2 + 3M) 



(2.8) 



Hence the phase-space variables can be written as A = 
(Aq, Ai, Ap), with the corresponding integration mea- 
sure as dt = d2(p+i)^_ 



C. Gaussian field operators 

The above results define a completely general Gaus- 
sian operator in terms of arbitrary bosonic annihilation 
and creation operators, without reference to the field 
involved. It is sometimes useful to compare this to a 
field theoretic notation, in which we explicitly use a 
coordinate-space integral to define the correlations. This 
provides a means to extend operator representation the- 
ory for fields 25, 26, 27, 28] to more general basis sets. In 
a quantization volume V , one can expand: 



*,(x) 
*!(x) 



-ik.x 



where the field commutators are 

"$,(x),$t (x')l ^%,5(x-x') 



(2.9) 



(2.10) 



With this notation, the quadratic term in the Gaussian 
exponent becomes: 



J2.n) 

where we have introduced the extended vector J[.(x) = 

(*, (^^''^ )i and 5$(x) = $(x) - ^(x) which is the op- 
erator ffuctuation relative to the coherent displacement 
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or classical mean field. If we index the extended vector 
as , where s = —1(1) for the first and second parts 
respectively, this Fourier transform can be written com- 
pactly as: 



1 



(2.12) 



The notation a ^ (x, y) indicates a functional matrix in- 
verse where: 



J g-\^,y)g{y,^')d'y = m^- , (2.13) 

and the relationship to the previous cross- variance matrix 
is that: 



1 



(2.14) 

In the standard terminology of many-body theory 
and field theory these field variances are generalized 
equal-time Green's functions, and can be written as: 



ct(x,x') = 



I(5(x, x') + n(x, x') m(x, x') 

m(x,x')+ K(x,x') + n^(x',x) 

(2.15) 

We shall show in the next section that these indeed cor- 
respond to field correlation functions in the case that the 
field state is able to be represented as a single Gaussian. 
More generally, one must consider a probability distribu- 
tion over different coherent fields and Green's functions 
or variances, in order to construct the overall density ma- 
trix. 



and using the eigenvalue property of coherent states, 
2 |z) = z |z), we find that: 



Tr 



A 



f^/d2*^zexp [-5z+qr^5z/2\ 



(3.2) 



The normalizing factor can now be recognized as the de- 
terminant expression arising in a classical Gaussian. For 
example, in the single-mode case, one obtains for the nor- 
malizing determinant that: 



(3.3) 



■\/(l + nY — mm+ 



We can thus calculate the value of the normalization from 
standard Gaussian integrals, as detailed in Appendix IbI 
provided a has eigenvalues with a positive real part. The 
result is: 



(3.4) 



Thus for A itself to correspond to a normalised density 
matrix, we must have = 1. In a general expansion of 
a density matrix, there may be terms which do not have 
this normalisation; with the proviso the average weight 
still be (il) = 1. This freedom of having different weights 
on different members of the ensemble provides a way of 
introducing gauge variables, which can be used to im- 
prove the efficiency of the stochastic sampling but which 
do not affect of the average result. The weight also allows 
calculations to be performed in which the trace of the 
density matrix is not preserved, as in canonical-ensemble 
calculations. 



III. GAUSSIAN EXPECTATION VALUES 

In order to use the Gaussian operator basis, a number 
of basic identities are needed. In this section, we derive 
relations between operator expectation values and mo- 
ments of the distribution. Such moments also show how 
the general Gaussian representation incorporates the pre- 
viously used methods. 



A. Gaussian trace 

The trace of a generalized Gaussian is needed to nor- 
malize the density matrix. The trace is most readily cal- 
culated by using a well-known coherent-state identity (^], 



Tr 

Here we define z 
tended vectors z 



A 



lAI 



-l2M„ 



rM 



(3.1) 



{zi,...,zm)- Next, introducing ex- 



(z,z*) 



*\T 



(z*,z), 5z 



B. Expectation values 

Given a density matrix expanded in Gaussian opera- 
tors, it is essential to be able to calculate operator expec- 
tation values. This can be achieved most readily if the 
operator O is written in antinormally ordered form, as: 



(3.5) 



Since the density matrix expansion is normally ordered 
by definition, the cyclic properties of a trace allows the 
expectation value of any antinormally ordered operator 
to be re-arranged as a completely normally ordered form. 
Hence, following a similar coherent-state expansion pro- 
cedure to that the previous subsection, we arrive at an 
expression analogous to the kernel trace, Eq H3.2() : 



O 



TrE,^,(at)p^,(a)] 
Tr[p] 
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/ p{x)o{\)ndx 
0{T 



(3.6) 



Here we have introduced an equivalence between the 
quantum expectation value (O), and the weighted prob- 



abilistic average (^0{ A)y . This is an antinormally or- 
dered c-number operator equivalence in phase space of 
0(A) ~ O, where the eigenvalue relations of coherent 
states are utilized to obtain: 



0(A) 



(oiz))- 



(3.7) 



Here (o(z))-r> represents the classical Gaussian average 
A 

of the c-number function o(z). In other words, all quan- 
tum averages are now obtained by a convolution of a 
classical Gaussian average with a width cta that depends 
on the kernel parameter A , together with a probabilis- 
tic average over A , with a width ap that depends on the 
phase-space distribution P( A ). The situation is depicted 
schematically in FigQ] 

Consider the first-order moment where O — a^^. This is 
straightforward, as o{z) — z^, and the Gaussian average 
of o(z) is simply the Gaussian mean a^: 



(3. 



More generally, to calculate the antinormally ordered 
moment o(a) — {a^j^a^2....a^^}, one must calculate the 
corresponding Gaussian moment o(z) — z^^z^^...z^^. 
This is most easily achieved by use of the moment- 
generating function for the Gaussian distribution in Eq. 
1)3. 7|) . which is 



XA{h A) 



t' a+t^g_t/2 



(3.9) 



where t = (ii, tM, i*, ^m) = (t,t*^)- General mo- 
ments of the Gaussian distribution are then given by: 



(3.10) 



t=o 



where it must be remembered that the adjoint vector 
t* is not independent of t. We note that averaging the 
moment-generating function over the distribution P{ X) 
gives the antinormal quantum characteristic function of 
the density operator: 



XA 



(t,t*) = Tr{pe**^e^**} 

P( A')r!e-'-+-=-/'dt . (3.11) 



This equation is an alternative way of (implicitly) defin- 
ing the Gaussian P distribution as a function whose gen- 
eralized Fourier transform is equal to the quantum char- 
acteristic function. 

As an example of a moment calculation, one obtains 
the c-number operator equivalence for general normally 
ordered quadratic term as: 



(3.12) 



where we have introduced the normally ordered covari- 
ance: 



= a-I. 



(3.13) 



Writing these out in more detail, we obtain the 
following central results for calculating normally ordered 
observables up to quadratic order: 





= ("i>p 












= {aiUj -1 


rriij ) p 
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= - 








+ m+)^ . (3.14) 



Comparing these equations with the schematic dia- 
gram in Fig n we see that, as expected from a convo- 
lution, the overall variance of any quantity is the sum of 
the variances of the two convolved distributions; that is, 
a = a\-\-ap. The results also support our interpretation 
given in section III Bl that n and m are, respectively, the 
normal and anomalous correlations that appear in many- 
body theory - except for the additional feature that we 
can now allow for distributions over these correlations. 
The expressions in the P-averages on the right-hand side 
are not complex-conjugate for Hermitian-conjugate op- 
erators, because the kernel A( A ) is generically not Her- 
mitian. Of course, after averaging over the entire distri- 
bution, one must recover an Hermitian density matrix, 
and hence the final expectation values of annihilation and 
creation operators will be complex-conjugate. Using the 
characteristic function, one can extend these to higher 
order moments via the standard Gaussian factorizations 
in which odd moments of fluctuations vanish, and even 
moments of fluctuations are expressed as the sum over 
all possible distinct pairwise correlations. 



C. Quantum field expectation values 

The results obtained above can be applied directly to 
obtaining the corresponding expectation values of nor- 
mally ordered field operators: 
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Representation 
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Wigner {W)\2] 
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Husimi (Q)f3l 
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Glauber-Sudarshan {P)\4] 
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Drummond-Gardiner (+P)[5] 
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a+ 











stochastic gauge [IQ., JJJ 
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Table I: Classification of commonly used single-mode operator 
representations in terms of parameters of the general Gaussian 
basis. 



*I(x) 

*.(x)$,(y) 

:*i(x)*J(y) : 
$l(x)$](y) 



(vl/,(x)) - (*.(x))p 

" ■ ^ = (*+(x))p 

= (*,(x)*j(y) +TO,j(x,y))p 
= (*.,(x)4'+(y)+n,,(x,y))^ 

= (vl'+(x)v|/+(y)+m+(x,y))^ . 

(3.15) 

These results show that in the field formulation of 
the Gaussian representation, the phase-space quantities 
riij (x, y) and (x, y) correspond to single-time Greens 
functions, analogous to those found in the propagator 
theory of quantum fields. 



D. Comparisons with other methods 

It is useful at this stage to compare these operator cor- 
respondences with the most commonly used previously 
known representations, as shown in Table For simplic- 
ity, this table only gives a single-mode comparison: 

In greater detail, we notice that 



• In the squeezed-state basis, the parameters n, m 
are not independent, as indicated in the table. The 
particle number n is a function n(|TO|)of the squeez- 
ing m. The exact relationship is given later. 

• The Gaussian family of representations is much 
larger than the traditional phase-space variety, be- 
cause we can allow other values of the cr^^ variance 
- for example, squeezed or thermal state bases. 
For thermal states, the variance corresponds to a 
Hermitian, positive-definite density matrix if Uij is 
Hermitian and positive definite, in which case riij 
behaves analogously to the Green's function in a 
bosonic field theory. In this case, a unitary trans- 
formation of the operators can always be used to 
diagonalize riij^ so that = UiSij. 

• For a general Gaussian basis, Gaussian operators 
that do not themselves satisfy density matrix re- 
quirements are permitted as part of the basis - 
provided the distribution has a finite width to com- 
pensate for this. This is precisely what happens, for 
example, with the well-known Q function, which 
always has a positive variance to compensate for 
the lack of fluctuations in the corresponding basis, 
which is Hermitian but not positive-definite. 

Distributions over the variance are also possible. It is 
the introduction of distributions over the variance that 
represents the most drastic change from the older distri- 
bution methods. It means that there many new operator 
correspondences to use. Thus, the covariance itself can 
be introduced as a dynamical variable in phase space, 
which can change and fluctuate with time. In this re- 
spect, the present methods have a similarity with the 
Kohn variational technique, which uses a density in co- 
ordinate space, and has been suggested in the context 
of BEC^24|. Related variational methods using squeezed 
states have also been utilized for BEG problems [23. By 
comparison, the present methods do not require either 
the local density approximation or variational approxi- 
mations. 



• If CT^i/ = S^f, these results correspond to the stan- 
dard ones for the normally ordered positive-P rep- 
resentation. 

• If we consider the Hermitian case of a* = 

as well, but with tr^i, = (n — l)(5^i/, where n ~ 
(s — l)/2, we obtain the 's-ordered' representation 
correspondences of Cahill and Glauber. 

• These include, as special cases, the normally or- 
dered Glauber-Sudarshan P representation (n = 
0), and the symmetrically ordered representation 
of Wigner {n = —1/2) 

• The antinormally ordered Husimi Q function is re- 
covered as the singular limit n —1 . 



IV. GAUSSIAN DIFFERENTIAL IDENTITIES 

An important application of phase-space representa- 
tions is to simulate canonical ensembles and quantum 
dynamics in a phase space. An essential step in this pro- 
cess is to map the master equation of a quantum density 
operator onto a Liouville equation for the probability dis- 
tribution P. The real or imaginary time evolution of a 
quantum system depends on the action of Hamiltonian 
operators on the density matrix. Thus it is useful to 
have identities that describe the action of any quadratic 
bosonic form as derivatives on elements of the Gaussian 
basis. These derivatives can, by integration by parts, be 
applied to the distribution P, provided boundary terms 
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vanish. The resultant Liouville equation for P is equiv- 
alent to the original master equation, given certain re- 
strictions on the radial growth of the distribution. When 
the Liouville equation has derivatives of only second or- 
der or less (and thus is in the form of a Fokker- Planck 
equation) , it is possible to obtain an equivalent stochastic 
differential equation which can be efficiently simulated. 

In general, there are many ways to obtain these iden- 
tities, but we are interested in identities which result in 
first-order derivatives, where possible. Just as for expec- 
tation values, this can be achieved most readily if the 
operator O is written in factorized form, as in Eq (|3.5|) . 

In this notation, normal ordering means: 



OA 



E 



(4.1) 



[a : Sa^A 



It follows that: 



an + (Jn 



A 



(4.4) 



(4.5) 



2. Quadratic products: 

Differentiating a determinant results in a transposed 
inverse, a result that follows from the standard cofactor 
expansion of determinants: 



d\g\ 



(4.6) 



We also need a notation for partial antinormal ordering: 
{d:A-^=Y,Ma):A:v,{a^), (4.2) 

i 

which indicates an operator product which antinormally 
orders all terms except the normal term ; A :. The Gaus- 
sian kernel A is always normally ordered, and hence we 
can omit the explicit normal-ordering notation, without 
ambiguity, for the kernel itself. 

In this section, for brevity, we use d /d A — 
(9/ri, i9/a, 9/q:+, 9/n, (9/m, (9/m+) to symbolize either 
d/dxi or —id/dyi for each of the i — 0, ..p complex 
variables A . This is possible since A( A ) is an ana- 
lytic function of A , and an explicit choice of deriva- 
tive can be made later. We first note a trivial identity, 
which is nevertheless useful in obtaining stochastic gauge 
equivalences between the different possible forms of time- 
evolution equations: 



(4.3) 



A. Normally ordered identities 

The normally ordered operator product identities can 
be calculated simply by taking a derivative of the Gaus- 
sian operator with respect to the amplitude and variance 
parameters. 



Linear products: 



Similarly, for the normalization factor that occurs in 
Gaussian operators: 



d\cj\ 



1 I I 

= -^^tiv 1^1 



(4.7) 



Hence, on differentiating with respect to the inverse 
covariance, we can obtain the following identity for any 
normalized Gaussian operator: 



-A 



: exp 



-So} -a-^ ■ Sa/2 



f_,„ ' Sa^Sal] A 



(4.8) 



Using the chain rule to transform the derivative, it fol- 
lows that a normally ordered quadratic product has the 
following identity: 



(5a^(5aj^A : = 



d 



dcTuu 



A 



A . (4.9) 



B. Antinormally ordered identities 

The antinormally ordered operator product identities 
are all obtained from the above results, on making use of 
the algebraic re-ordering results in the Appendix: 



The result for linear operator products follows directly 
from differentiation with respect to the coherent ampli- 
tude, noting that each amplitude appears twice in the 
exponent: 



d 



A = 



d n 



exp 



-(5aV"^(5a/2 



1. Antinormal linear products 

Antinormally ordered linear products can be trans- 
formed directly to normally ordered products. Hence, 
from the Appendix and Eq H4.5(l . we obtain: 



: A :| = : [a,, - a^lSa^] A : 
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d 



d 



A 



A, 



(4.10) 



where we recall from Eq H3.13|l that the normaUy ordered 
covariance is defined by: = ^ — /. 



2. Quadratic products with one antinormal operator 

This calculation follows a similar pattern to the previ- 
ous one: 



|(5a^ : 5alA :| 



5a,, 



d 
dal 

[St_cu + {St_cp - (J-^)SalSap] A 
d 



A . (4.11) 



3. Quadratic products with two antinormal operators 

We first expand this as the iterated result of two re- 
orderings, then apply the result for a linear antinormal 
product to the innermost bracket: 

= {Sa^ [5.p-<7;.'] : <5StA :}(4.12) 

Next, the result above for one antinormal operator is 
used: 



[Sa^Sal : A :| = [S^p - a^^] 

X (Tpp + 2{apa — 6pa)<7t3p 

d 



d 



A 



da-pa 



da/3a_ 

A . (4.13) 



C. Identities in matrix form 

The different possible quadratic orderings can be writ- 
ten in matrix form as 



: aa'A 



[a ■■a' A :| 



: aa^A : Aaa^ 

at^atA at^Aa^ 

aAat aa^A 

Aat^^at {at^^Aa^} 

aa^A aAa^ 

at^Xat |at^:a^A:} 



.(4.14) 



With this notation, all of the operator identities can 
be written in a compact matrix form. The resulting set 
of differential identities can be used to map any possible 
linear or quadratic operator acting on the kernel A into 
a first-order differential operator acting on the kernel. 

For this reason, the following identities are the central 
result of this paper: 



A 

: aA : 

{aA} 

: 6a^5a^A : 
|,5a:(5a^A:} 
\^SaSa<A^ 



= aA 



dA 



aA + a 



N 



OA 
da+' 



T o 9A 
a A + 2g_—a 

— —OCT — 



dA 



g-A + 2q!' —a , 

— — dg_— 

^N'l: , r,N ^A « 
g_ iv -|- Zg_ — — g_ 

— — da — 



(4.15) 



Here the derivatives are defined as 



d_ 

da 

d 



da^ 




d_ / d 
da. ' \9q;+ 

d/da+ \ 
{d/dccf ) 

d 



da: 



(4.16) 



p,u 



up 



It should be noted that the matrix and vector derivatives 
involve taking the transpose. We note here that for no- 
tational convenience, the derivatives with respect to the 
ap,i, are formal derivatives, calculated as if each of the 
(^v^p were independent of the others. With a symmetry 
constraint, the actual derivatives of A with respect to 
any elements of n or any off-diagonal elements of m will 
differ from the formal derivatives by a factor of two. For- 
tunately, because of the summation over all derivatives 
in the final Fokker-Planck equation, the final results are 
the same, regardless of whether or not the symmetry of 
ap^n is explicitly taken into account at this stage. 

The quadratic terms can also be written in a form with- 
out the coherent offset terms in the operator products. 
This is often useful, since while the original Hamiltonian 
or master-equation may not have an explicit coherent 
term, terms like this can arise dynamically. The following 
result is obtained: 



: aa^A 



dA 
da — 



dA 

r 

-da' 



-a 
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(aa + cr) A + 2a— a, 
^ —' —aa — 



(a :a'''A:l = o^q_^ 
y J da- — 



dA 



da— — da^ 



OA 



(32+ + a)A + 2ct^— CT, 
^ —' — da — 



<aa'A} = a— a +a ^r—qP 
L J oa— — oa+ 



.OA 



(aa+ + a'^)A + 2a^^a^ (4.17) 
^ — ' — oa — 



One consequence of these identities is that the time 
evolution resulting from a quadratic Hamiltonian can 
always be expressed as a simple first-order differential 
equation, which therefore corresponds to a determinis- 
tic trajectory. This relationship will be explored in later 
sections: it is quite different to the result of a path in- 
tegral, which gives a sum over many fluctuating paths 
for a quadratic Hamiltonian. Similarly, the time evolu- 
tion for cubic and quartic Hamiltonians can always be 
expressed as a second-order differential equation, which 
corresponds to a stochastic trajectory. 



D. Identities for quantum field operators 

The operator mappings can also be succinctly written 
in the field-theoretic notation as 



derivatives have been defined as 



d _ 1 



isk.x 



da 



d 



9*+(x) 



d 



, — i(s'k-x— sk'-x') 



dajsj's'ix,^') ^ k daw'jsMj's'' 
Again we have the convention for matrix derivatives that 
f d \ d 



da{-x., x') 



js,j's' 



dai 



.(x,x') 



V. EXAMPLES OF GAUSSIAN OPERATORS 

This section focuses on specific examples of Gaussian 
operators, and relates them to physically useful pure 
states or density matrices. We begin by defining the class 
of Gaussian operators that correspond to physical density 
matrices, before looking at examples of specific types of 
states that can be represented, such as coherent, squeezed 
and thermal. In each of these specific cases, the conven- 
tional parametrization can be analytically continued to 
describe a non-Hermitian basis for a positive representa- 
tion. We show how these bases include and extend those 
of previously defined representations, and calculate the 
normalization rules and identities that apply in the sim- 
pler cases. 



:\E:.(x)A: = £(x)A + /dVCT(x,x') 



dA 



{$(x)a} = £(x)A + y"dV^^(x,x' 



a*.+ (x') 
dA 



(5l(x)5$(x')^A : = o:(x,x')A + 2 / /dV'rf 



X 



o:(x,x")- 



dA 



I III i\ 
:o:(x ,x) 



{5i(x) : (5$(x')^A :} 



|ji(x)(5$(x')tA} 



a^(x",x"') = 
o:(x,x')A + 2 / ld\"d\"' X 



a^(x,x") 



dA 



r^(x'",x') 



a^(x",x"')= 

CT^(x,x')A + 2 lld?x"d\"' X 



^^(x,x") 



dA 



A. Gaussian density matrices 

A Gaussian operator can itself correspond to a physi- 
cal density matrix, in which case the corresponding dis- 
tribution is a delta function. This is the simplest pos- 
sible representation of a physical state. Gaussian states 
or physical density matrices are required to satisfy the 
usual constraints necessary for any density matrix: they 
must be Hermitian and positive definite. From the mo- 
ment results of Eq (|3.14|l , the requirement of Hermiticity 
generates the following immediate restrictions on the dis- 
placement and covariance parameters: 



N I III i\ 

a (x ,x). 



ao:(x",x"')^ 



a. 
n 



(5.1) 



(4.18) 



where the vector quantum fields and covariances are as 
defined in section Hi CI The normal field correlation ma- 
trix is a;^(x, x") — ct(x, x") — /(5(x, x') and the functional 



In addition, there are requirements due to positive- 
definiteness. To understand these, we first note that 
when n is Hermitian, as it must be for a density matrix, 
it is diagonalizable via a unitary transformation on the 
mode operators. Therefore, with no loss of generality, we 
can consider the case of diagonal n, i.e., n^j =nkSkj- The 
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positive-definiteness of the number operator then means 
that the number eigenvalues are real and non-negative: 



rik > 



(5.2) 



In the diagonal thermal density matrix case, but with 
squeezed correlations as well, satisfying the density 
matrix requirements means that there are additional 
restrictions^!. Consideration of the positivity of prod- 
ucts like XkjXlj where: = fia^ + means that one 

must also satisfy the inequalities nfc(l -I- nj) > \mkj\'^- 
This implies a necessary lower bound on the photon num- 
ber in each mode: 



rik > n{\mkk\) = Vlm/cfep + 1/4-1/2 



(5.3) 



Examples of Gaussians of this type are readily obtained 
by first generating a thermal density matrix, then apply- 
ing unitary squeezing and/or coherent displacement op- 
erations, which preserve the positive definite nature of 
the original thermal state. This produces a pure state if 
and only if the starting point is a zero-temperature ther- 
mal state or vacuum state. Hence, the general physical 
density matrix can be written in factorized form as: 



Here 



Ati^(n) = : exp -2+ [1 + n] ^ a 



(5.4) 



(5.5) 



is a thermal density matrix completely characterized by 



Tr 



: Ath{n) 



where 



its number expectation: n 

n must be Hermitian for the operator to correspond to a 
physical density matrix. We show the equivalence of this 
expression to the more standard canonical Bose-Einstein 
form in the next section. 

The unitary displacement and squeezing operators are 
as usually defined in the literature: 



and: 



(5.6) 



(5.7) 



where the vector a is, as before, the coherent displace- 
ments for each mode. The symmetric matrix ^ gives the 
angle and degree of squeezing for each mode, as well as 
the squeezing correlations between each pair of modes. 

In table iQTJ, we give a comparison of the Gaussian 
parameters found in the usual classifications of physical 
density matrices of bosons, for a single-mode case. 



B. Thermal operators 

1. Physical states 



Physical state 


Q, 


a 


Q + 


n 


m 


m+ 


Vacuum state 


1 

















Coherent state 


1 


a 


a* 











Thermal 


1 








n > 








Squeezed vacuum 


1 








n{\m\) 


m 


m* 


Squeezed coherent 


1 


a 


a* 


n{\m\) 


m 


* 

m 


Squeezed thermal 


1 


a 


a* 


n > n(\m\) 


m 


m* 



Table II: Parameters of single-mode Gaussian density matri- 
ces of bosons. 



cal form aslSlI: 



(5.8) 



where (j)^ — Sk/kT . Here the modes are chosen, with 
no loss of generality, to diagonalize the free Hamiltonian 
with mode energies , and for the case of massive bosons 
we have included the chemical potential in the definition 
of the energy origin. To show how this form is related 
to the normally ordered thermal Gaussian A-Q^{ri) of Eq 
(|5.5|) , we simply note that since n is hermitian, it can be 
diagonalized by a unitary transformation. The resulting 
diagonal form in either expression is therefore diagonal 
in a number state basis and is uniquely defined by its 
number state expectation value. 

Clearly, one has for the usual canonical density matrix 
that: 



\Pth\ 



exp I 



5fef^fc 



(5.9) 



while it is straight-forward to show that the correspond- 
ing normally ordered expression is a binomial: 



(n|Ath(n)|n)=n[l 



- rik] ^ 



1 



1 



1 + nfe 



(5.10) 

As one would expect, these expressions are identical pro- 
vided one chooses the standard Bose-Einstein result for 
the thermal occupation as: 



nk 



1 



(5.11) 



These results also show that when n = one has a 
vacuum state, corresponding to a bosonic ground state 
at zero temperature. In summary, the normally ordered 
thermal Gaussian state is completely equivalent to the 
usual canonical form. 



2. Generalized thermal operators 



It is conventional to write the bosonic thermal density 
operator for a non-interacting Bose gas in grand canoni- 



A simple non-Hermitian extension of the thermal 
states can be defined as an analytic continuation of the 
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usual Bose-Einstein density matrix for bosons in ther- 
mal equilibrium. We define a normally ordered thermal 
Gaussian operator as having zero mean displacement and 
zero second or fourth quadrant variance: 



A(r2,0,0,n,0,0) 



exp 



(I- 



(5.12) 

Such operators are an analytic continuation of previously 
defined thermal bases, and are related to the thermo- field 
methods [s^. 

As well as the usual Bose-Einstein thermal distribu- 
tion, the extended thermal basis can represent a variety 
of other physical states. As an example, consider the gen- 
eral matrix elements of an analytically continued single- 
mode thermal Gaussian operator in a number-state basis, 
with 1 + 1/n = exp [(f)] = exp [(r -I- itp)]. These are: 



(n|Ath(^)l"'> = ("I [1 - e-^] exp [-H l*^') 

= Kn' [l - e^"*] exp [-n{r + iif))] . 

(5.13) 

Now consider the following single-mode density matrix: 



1 /■ _ -1 _ 

p = — / [1 - e""^] exp [no (r i V)] {n)d'ip . 
^'^ Jo 

(5.14) 

Taking matrix elements in a number-state basis gives: 



{n\p\n) = 



2ir 



'{r+i^){n-no) 



dip 



(5.15) 



This effectively Fourier transforms the thermal opera- 
tor on a circle of radius |no| around the origin, thereby 
generating a pure number state with boson number equal 
to riQ. Thus, extended thermal bases of this type are 
certainly able to represent non-Gaussian states like pure 
number states. Nevertheless, they cannot represent co- 
herences between states of different total boson number. 



3. Thermal operator identities 

The operator identities for the thermal operators are a 
subset of the ones obtained previously. There are no use- 
ful identities that map single operators into a differential 
form, nor are there any for products like aicij. However, 
all quadratic products that involve both annihilation and 
creation operators have operator identities. 

With this notation, and taking into account the fact 
that differentiation with respect to n now explicitly pre- 
serves the skew symmetry of the generalized variance, the 
operator identities can be written: 



A 

: aa^A : 
{a : a^A :} 
{:aA : at} 

{aa^A} 



4^ 



dA 



(l + n)A+(l + n) — (1 + n) 

1 + n A + n— 1 + n , 
on 

/-, ^dA 
(l + n)A + (l + n)— n , 

aA 

nA + n— — n . 
an 



(5.16) 



C. Coherent projectors 

1. Physical states 

Next, we can include coherent displacements of a ther- 
mal Gaussian in the operator basis. This allows us to 
compare the Gaussian representation with earlier meth- 
ods using the simplest type of pure-state basis, which is 
the set of coherent states. These have the property that 
the variance in position and momentum is fixed, and al- 
ways set to the minimal uncertainty values that occur in 
the ground state of a harmonic oscillator. 

In general, we consider an Af-mode bosonic field. In an 
M-mode bosonic Hilbert space, the normalized coherent 
states I a) are the eigenstates |a) of annihilation opera- 
tors a with eigenvalues a. The corresponding Gaussian 
density matrices are the coherent pure-state projectors: 



Ada) 



(5.17) 



which are the basis of the Glauber-Sudarshan P- 
representation. To compare this with the Gaussian no- 
tation, we re-write the projector using displacement op- 
erators, as: 



Ac (a) = e°^-° |0) (0|e°*-°-l" 



(5.18) 



Since the vacuum state is an example of a thermal Gaus- 
sian, and the other terms are all normally ordered by con- 
struction, this is exactly the same as the Gaussian opera- 
tor A(l, a, a*, 0, 0, 0) . In other words, if we restrict the 
Gaussian representation to this particular subspace, it is 
identical to the Glauber-Sudarshan P-representation 4| . 
This pioneering technique was very useful in laser physics, 
as it directly corresponds to easily measured normally or- 
dered products. It has the drawback that it is not a com- 
plete basis, unless the set of distributions is allowed to in- 
clude generalized functions that are not positive-definite. 

Other examples of physical states of this type are 
the displaced thermal density operators. These physi- 
cally correspond to an ideal coherently generated bosonic 
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mode from a laser or atom laser source, together with a 
thermal background. They can be written as: 

Ac{a,n) — A(l, a, a*, n, 0, 0) 

= e°* """°*Atl^(ri:)e°' °^""* ° . (5.19) 

2. Generalized coherent projectors 

There are two ways to generalize the coherent projec- 
tors into operators that are not density matrices: either 
by altering the thermal boson number n so it does not 
correspond to a physical state, or by changing the dis- 
placements so they are not complex conjugate to each 
other. 

The first procedure is the most time-honored one, since 
it is the route by which one can generate the classical 
phase-space representations that correspond to different 
operator orderings. The Wigner ^, Q-function's^, and 
s-ordered|H bases are very similar to Gaussian density 
matrices, except with negative mean boson numbers: 

Aw (a) A(1,q;,q;*, -//2,0,0) 
Aqia) = A(l,a,a*,-/,0,0) 
A, (a) = A(l,a,a*,7(s- l)/2,0,0) . (5.20) 

As pointed out in the previous subsection, it is also pos- 
sible to choose n to be non-Hermitian, which would al- 
low one to obtain representations of coherently displaced 
number states. However, there is a problem with this 
class of non-normally ordered representations. Gener- 
ically, they have a restricted set of operator identities 
available, and typically lead to Fokker-Planck equations 
of higher than second order - with no stochastic equiva- 
lents - when employed to treat nonlinear Liouville equa- 
tions. 

Another widely used complete basis is the scaled 
coherent-state projection operator used in the positive-P 
representation 9|| and its stochastic gauge extensions p^: 

AAn,a,^)^n\^. (5.21) 

Here we have introduced (3* as a vector amplitude for 
the coherent state |/3*), in a similar notation to that used 
previously. 

This expansion has a complex amplitude f2, and a dy- 
namical phase space which is of twice the usual classical 
dimension. The extra dimensions are necessary if we wish 
to include superpositions of coherent states, which give 
rise to off-diagonal matrix elements in a coherent state 
expansion. To compare this with the Gaussian notation, 
the projector is re-written using displacement operators, 
as: 

Apin,a,l3) = Sle°^'^-« |0) (0|e^-°-'3-° 

= A(17,a,/3,0,0,0). (5.22) 



This follows since the vacuum state is an example of 
a thermal Gaussian, and the other terms are all nor- 
mally ordered by construction. From earlier workQ, it 
is known that any Hermitian density matrix p can be 
expanded with positive probability in the over-complete 
basis Ap, and it follows that the same is true for A( A ). 



The effects of the annihilation and creation operators 
on the projectors are obtained using the results for the 
actions of operators on the coherent states, giving: 



A = 




A 






aA = 


aA 








S^A = 




/3 + 


d ' 
da 


A 




Aa ^ 




a + 


d ' 


A 




Aat = 


(3 A . 






(5.23) 



Note that here one has a — 0, and thus all 
the antinormally ordered identities have just coherent 
amplitudes without derivatives, in agreement with the 
general identities obtained in the previous section. In 
treating nonlinear time-evolution, this has the advantage 
that some fourth-order nonlinear Hamiltonian evolution 
can be treated with only second-order derivatives, which 
means that stochastic equations can be used. In a similar 
way, one can treat some (but not all) quadratic Hamilto- 
nians using deterministic evolution only. The fact that all 
derivatives are analytic - which is possible since (3 ^ a* 
- is an essential feature in obtaining stochastic equations 
for these general cases 0. 

D. Squeezing projectors 

1. Physical states 

The zero-temperature subset of the Gaussian den- 
sity operators describe the set of minimum uncertainty 
states, which in quantum optics are the familiar squeezed 
states 33, 34]. These are most commonly defined as the 
result of a squeezing operator on a vacuum state, followed 
by a coherent displacement: 

Asq(a, 1, 0) = D{a)S{i} |0) (0| S{~i)D{~a) . (5.24) 

The action of the multi-mode squeezing operator on an- 
nihilation and creation operators is to produce 'anti- 
squeezed operators: 

b = 5(1)25^(1) = fia + i^a^^ 

fot^ = 5(|)at^^t(0 - M*at^ + u*a , (5.25) 

where the Hermitian matrix /x(^) and the symmetric ma- 
trix i''(^) are defined as multi-mode generalizations of hy- 
perbolic functions [35I l3^ : 
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= cosh(lll) 
_ sinh(l^l) 



(5.26) 



Note that and v obey the hyperbohc relation /x/x — 
i/i/* = I and have the symmetry property ij,~^v — 
=i'*ij,~^. In the physics of Bose-Einstein con- 
densates, b and bt are just the Bogohubov annihilation 
and creation operators for quasi-particle excitations. 

The Bogoliubov parameters provide a convenient way 
of characterizing the minimum uncertainty Gaussian op- 
erators. We therefore need to relate them to the param- 
eters in the Gaussian covariance matrix. First consider 
the antinormal density moment for a squeezed state: 

{ad)) = Tr{aStAsq(a,tO)} 

= (0| S{-$,)D{-a)Ud)D\-OL)S\-£,) |0) 
= (0| [na - jyflt + a) (a^ fi - av* + a*) |0) 
= OLCx* + /x/x. (5.27) 

Similarly, the anomalous moments are: 



(5.28) 



Comparing these moments to those of the general Gaus- 
sian state [Eq. (|3.14|) ]. we see that 



n = uu 
m = fjLi/ 
m = u fi. 



(5.29) 



The relationship between the different parameteriza- 
tions can be written in a compact form if we make the 
definitions 




(5.30) 



in terms of which the relations are 

1 2 1. 
2= 2- 



fi = exp 



(5.31) 



One implication of this relation is that, just as /x is not in- 
dependent of I/, so too n is not independent of m for the 
squeezed state. From the hyperbolic relation, we see that 



= m* (1 + 'n)~^ m. The determinant of the covari- 
ance matrix, required for correct normalization, reduces 
to the simpler form: 



a = |l + n| = |/x|- 



(5.32) 



This set of diagonal squeezing projectors forms the ba- 
sis that has previously been used to define squeezed-state 
based representations'^. Because the basis elements in 
such bases are not analytic and the resultant distribution 
not always positive, these previous representations suffer 
from the same deficiency as the Glauber-Sudarshan P 
representation, (as opposed to the positive-P representa- 
tion), i.e. the evolving quantum state can not always be 
sampled by stochastic methods. 



2. Generalized squeezing operators 

A non-Hcrmitian extension of the squeezed-state basis 
[Eq. (|5.24|) ] can be formed by analytic continuation of its 
parameters, i.e. by a replacement of the complex con- 
jugates of a and ^ by independent matrices: cx* 
and ^* — > In the Bogoliubov parametrization, this 
is equivalent to the replacement u* f+and to /x be- 
ing no longer Hermitian. These non-Hermitian opera- 
tors are in the form of off-diagonal squeezing projectors 
and constitute the basis of a positive-definite squeezed- 
state representation. They include as a special case 
{i/ — i/+ — 0, fi — I) the kernel of the coherent-state 
positive-P expansion. Thus the completeness of the more 
general representation is guaranteed by the completeness 
of the coherent-state subset, and we can always find a 
positive-P function for any density operator by using the 
coherent-state based representation. 



E. Thermal squeezing operators 

Mixed (or classical) squeezed states are generated by 
applying the squeezing operators to the thermal kernel, 
rather than to the vacuum projector: 



Asq(0,en)-5(OAth(n)^^(-O 



(5.33) 



In this way, a pure or mixed Gaussian state of arbitrary 
spread can be generated. 

Once again, we can relate the covariance parameters 
characterizing the final state to the thermal and squeez- 
ing parameters by comparing the moments: 

(aat) = Tr{aatA,q(0,£,n)} 

= Tr |(/xa - iva^) (aV - az/*) A^^(n)| 
= n {n + I) fi + vTi'^ u* , (5.34) 
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since there arc no anomalous fluctuations in a thermal 
state. Similarly, the squeezing moments are 

<|22"^) = —fi {n + I) V ~ unF (1* 
(a^^a^) = -ij,*n'^u* -u* (rT+I)fM. (5.35) 

Thus the two parameterizations are related by 

n = finfi + i> (w^ + /) I/* 
m = —fi {n + I)v — uriF IX* 
m* = -n*rfv>*-u*{n + I)n, (5.36) 

which can be written in a compact form as 



£ = 11 + ^ij t + 
where the thermal matrix is defined as 



n 



n 




(5.37) 



(5.38) 



of the density matrix. To understand why it is useful to 
treat both types of evolution with the same representa- 
tion, we recall that the quantum theory of experimen- 
tal observations generally requires three phases: state 
preparation, dynamical evolution, and measurement. It 
is clearly advantageous to carry out all three parts of the 
calculation in the same representation, in order that the 
computed trajectories and probabilities are compatible 
throughout. Many-body state preparation is nontrivial, 
and often involves coupling to a reservoir, which may re- 
sult in a canonical ensemble. This can be computed using 
imaginary time evolution, as explained below. Dynamical 
evolution typically requires a real-time master equation, 
while the results of a measurement process are operator 
expectation values, which were treated in Section (III). 



A. Operator Liouville Equations 

Either real or imaginary time evolution occurs via a 
Liouville equation of generic form: 



As in the cases for the other bases, these squeezed ther- 
mal states can be analytically continued to form a non- 
Hermitian basis for a positive-definite representation. 
Such a representation would be suited to Bose-condensed 
systems, which have a finite-temperature (thermal) char- 
acter as well as a quantum (squeezed, or Bogoliubov) 
character. Furthermore, the lack of a coherent displace- 
ment is natural in atomic systems, where superpositions 
of total number are unphysical. 



F. Displaced thermal squeezing operators 

Finally, the most general Gaussian density matrix is 
obtained as stated earlier, by coherent displacement of a 

squeezed thermal state: 

A,q(a,C,n) =5(a),S(0Ath(n)|5(-4)5(-a) . (5.39) 

In this way, a pure or mixed Gaussian state of arbitrary 
location as well as spread can be generated. In terms 
of the normally ordered Gaussian notation, the displace- 
ment and covariance of this case are given by: 



o; = /Lt I n 




« = J • (5-40) 



VI. TIME EVOLUTION 

The utility of the Gaussian representation arises when 
it is used to calculate real or imaginary time evolution 



-m = L{p{t)) , 



(6.1) 



where the Liouville super-operator typically involves pre- 
and post-multiplication of p by annihilation and creation 
operators. There are many examples of this type of equa- 
tion in physics (and, indeed, elsewhere). We will con- 
sider three generic types of equation here: imaginary time 
equations used to construct canonical ensembles, unitary 
evolution equations in real time, and general non-unitary 
equations used to evolve open systems that are coupled 
to reservoirs. 

Wc often assume that initially, the steady-state density 
matrix is in a canonical or grand canonical ensemble, of 
form: 



-TH/h 



(6.2) 



where Pu{t) is un- normalized, r = h/kT, and we can 
include any chemical potential in the Hamiltonian with- 
out loss of generality. If this is not known exactly, the 
ensemble can always be calculated through an evolution 
equation in r, whose intial condition is a known high- 
temperature ensemble. This equation can also be ex- 
pressed as a master equation, though not in Lindblad 
form. The resulting equation in 'imaginary time', or r, 
can be written using an anti-commutator: 



H,p^ 



(6.3) 



Here the initial condition is just the unit operator. 

By comparison, the equation for purely unitary time 
evolution under a Hamiltonian H is: 



H,p 



(6.4) 
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More generally, one can describe either the equilibra- 
tion of an ensemble, or non-equilibrium behavior via a 
master equation representing the real-time dynamics of 
a physical system. Equations for damping via coupling 
of a system to its environment must satisfy restrictions 
to ensure that p remains positive-definite. In the Marko- 
vian limit, the resulting form is known as the Lindblad 
form 37]: 



dp 
dt 



K 



(6.5) 



which consists of a commutator term involving the Her- 
mitian Hamiltonian operator iJ, as well as damping 
terms involving an anti-commutator of the arbitrary op- 
erators Ok- 



B. Phase-space mappings 

While the general operator equations become exponen- 
tially complex for large numbers of particles and modes, 
the use of phase-space mappings provides a useful tool 
for mapping these quantum equations of motion into a 
form that can be treated numerically, via random sam- 
pling techniques. 

Using the operator identities in Eq. I|4.15() . one can 
transform the operator equations in any of these three 
cases into an integro-differential equation: 



dt 



(FX 



(6.6) 



where the differential operator Ca is of the general form 



Ca^U + A,dj + -D,jd^dj, 



(6.7) 



with derivative operators to the right, and i, j = 0, 
for the case of a p— parameter Gaussian. We only con- 
sider cases where terms with derivatives of order higher 
than two do not appear, which implies a restriction on 
the nonlinear Hamiltonian structure. 

We next apply partial integration to Eq. (|6.6|l . which, 
provided boundary terms vanish, leads to a Fokker- 
Planck equation for the distribution: 



d_ 

dt 



PiX,t) = CNP{X,t), 



(6.8) 



where the differential operator Cn has derivatives to the 
left: 

Cn = U~d,Aj + ^d,djD,j. (6.9) 

Such Fokker-Planck equations have equivalent path- 
integral and stochastic forms, which can be treated with 
random sampling methods. 



For example, in the Hamiltonian case, if the original 
Hamiltonian iJ7v(S, a^) is normally ordered (annihilation 
operators to the right) , then for a positive-P representa- 
tion one can immediately obtain: 



-N 



1 

ih 



[HNia,l3~do.)-HNif3,a^dfi)] . (6.10) 



With the use of additional identities in fl to eliminate 
the potential term U, the Fokker-Planck equation can be 
sampled by stochastic Langevin equations for the phase- 
space variables. Note that this potential term only arises 
with imaginary-time evolution. The first-order deriva- 
tive (drift) terms in the Fokker-Planck equation map to 
deterministic terms in the Langevin equations, and the 
second-order derivative (diffusion) terms map to stochas- 
tic terms. To obtain stochastic equations, we follow the 
general stochastic gauge technique 10], which in turn is 
based on the positive-P method. 

To simplify notation, we have left the precise form 
of the derivatives in the Fokker-Planck equation as yet 
unspecified. Different choices are possible because the 
Gaussian operator kernel is an analytic function of its pa- 
rameters. The standard choice in the positive-P method, 
obtained through the dimension-doubling technique'^, is 
such that when the equation is written in terms of real 
and imaginary derivatives, all the coefhcients are real 
and the diffusion is positive-definite. This ensures that 
stochastic sampling is always possible. Other choices are 
also possible and useful if analytic solutions are desired. 

The structure of the noise terms in the stochastic equa- 
tions is given by the noise-matrix B, which is defined as 
a p X p' complex matrix square root: 



D = BB' 



(6.11) 



Since this is non-unique, one can introduce diffusion 
gauges from a set of matrix transformations U[f(A)] 
with UU"^ = I. It is also possible to introduce arbi- 
trary drift gauge terms g which are used to stabilize the 
resulting stochastic equations: 



dfi 

'dt 
dXj 

dt 



A, + Py [0(0- 5j], (^,J>0) . (6.12) 



These are Ito stochastic equations with noise terms de- 
fined by the correlations: 



m)QAi'))^5{t-t')5, 



(6.13) 



We note here that the use of stochastic equation sam- 
pling as described here represents only one possible way 
to sample the underlying Fokker-Planck equation. Other 
ways are possible, including the usual Metropolis and 
diffusion Monte-Garlo methods found in imaginary-time 
many-body theory. 
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In the remainder of this section, we consider quadratic 
Hamiltonians, or master equations. We show that under 
the Gaussian representation, these give rise to purely de- 
terministic or 'drift' evolution. We first treat the thermal 
case, then derive an analytic solution to the dynamics 
governed by a general master equation that is quadratic 
in annihilation and creation operators. Following this are 
several examples which show how the analytic solution 
can be applied to physical problems. While these exam- 
ples can all be treated in other ways, they demonstrate 
the technique, which will be extended to higher-order 
problems subsequently. 



1. Imaginary time evolution 

We consider this case in detail, even though it is rel- 
atively straightforward, because it gives an example of 
phase-space evolution which would require diffusive or 
stochastic equations using previous methods. The equa- 
tion in 'imaginary time', or r can be written using an 
anti-commutator. Since we are only considering linear 
evolution here, the relevant Hamiltonian is always diag- 
onalizable, and can be written as: 



(6.15) 



C. General quadratic master equations 

Any quadratic master equation can be treated exactly 
with the Gaussian distribution. To demonstrate this, we 
can cast any quadratic master equation into the form 



+Buf, {afj,al -.p:} 



+ {a^ : alp:] , 
= -h ^ : ap : +3^^^ {a:p:] 

g}p: +M^^g^ : p :| -f-£|a : a^yO :| 

(6.14) 



Tr 



A : a 



where the trace is a matrix structural operation (indi- 
cated by the double underline), not a trace over the op- 
erators. Here A^^'' is a real number, while ^'^'and -B'^' 
are complex column vectors with the generalised Hermi- 
tian property of A^^>^ = A^'^+ , B^^>^ = B^^^+ . 

The quadratic terms A, B_ and C_ are complex-number 
matrices that have the implicit superscript (2) dropped 
for notational simplicity. By construction, A and B_ pos- 
sess all the skew symmetries of a;: A — and B = B^ 
, i.e. they are Hermitian in the generalized sense defined 
earlier. The matrix C_ possesses only some of these skew 
symmetries, namely that the upper right and lower left 
blocks are each symmetric. Furthermore, the Hermitic- 
ity of the density operator requires that the matrix Her- 
mitian conjugate is equal to the generalized Hermitian 
conjugate: ^'^ = ^+ and = £+. 

By expanding p in the general Gaussian basis and ap- 
plying the operator identities in Eq. H4.15|l . we obtain a 
Liouville equation for the phase-space distribution P that 
contains only zeroth- and first-order derivatives. Since 
this can be treated by the method of characteristics, the 
time-evolution is deterministic: every initial value cor- 
responds uniquely to a final value, without diffusion or 
stochastic behavior. This can also be solved analytically, 
since the time-evolution resulting from a quadratic mas- 
ter equation is linear in the Gaussian parameters A . 



Next, we need to cast the un-normalized density oper- 
ator equation: 



H,Pu 



(6.16) 



into differential form. All the terms are of mixed form, 
including both normal and antinormal ordered parts, so 
the master equation can be written as: 



d ^ 



g{a :a^Pu:}]+A^°^Pu , (6.17) 



where: 



u) 

w'^ 

-Trw . 



(6.18) 



Using the identities in Eq H4.15|l . one finds that the 
corresponding differential operator to be 



(6.19) 



This leads to the following equation for the distribution, 
after integration by parts (which requires mild restric- 
tions on the initial distribution): 



dP V- 

a7 = 



d 
duk 



nk 



(1 + rife) P .(6.20) 



Solving first-order Fokker-Planck like equations in this 
form leads to deterministic characteristic equations: 



k 

hk = -ujknk (1 + rik) ■ 



(6.21) 



Integrating the deterministic equation for the mode oc- 
cupation rifc leads to the Bose-Einstein distribution also 
encountered in Eq (|5.11|l : 



nk 



1 



(6.22) 
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The weighting term occurs because this method of ob- 
taining a thermal density matrix resuhs in an un- 
normaUzed density matrix with trace equal to ^{t). 
From integration of the above equation one finds, as ex- 
pected from Eq H5.8|l . that: 



(6.23) 



where the steady state with coherent transients is given 
by: 



(6.28) 



For a quadratic master equation in Lindblad form, the 
Hamiltonian and damping operators can be expressed as 



2. Real time evolution 

In the Lindblad form of a master equation which is 
relevant to real time evolution, further restrictions apply 
to its structure than just the symmetries given above. 

The preservation of the trace of p in real-time master 
equations requires that A^^) = -B'-^\ In addition, we 
require that = {A + = and that the 

matrix sum D_=A + B_+C_ is anti-skew-symmetric : 
= —D^. The resultant differential equation for P is 
simplified by the fact that most of the symmetric terms 
from the identities are multiplied by the antisymmetric 
D , and thus give a trace of zero. In particular, the zeroth- 
order terms will sum to zero, leaving 



dP 
'dt 



oa 



P, 



(6.24) 

where ^=2A + Q= -^R - £+ . 

A Fokker-Planck equation such as this that contains 
only first-order derivatives describes a drift of the distri- 
bution and can be converted into equivalent deterministic 
equations for the phase-space variables: 



H 



Ok - 



Tr H : aa^ 



Q-kQ. 



(6.29) 



where, in block form. 



H 



Ok = 



n* — 

^K ~ 



O 



(1) 

K 

(2) 



^K 



O 



(1), 
K 



a 



(2)* 
K 



(6.30) 



Thus the coefficients of density terms : aa^ : appear in 
J?'^^ and the coefficients of squeezing terms aa^ appear 
in H^^\ The commutator term and each of the damping 
terms will provide a contribution to the matrices A, _B, 
and C_, which we can label respectively as A^, A^, etc. 
The contributions from the Hamiltonian term are thus: 



Ah = 



C 



H 



\ 

-iH(^> ) ' 

-2iff(i) 

2ii?(i)^ 



(6.31) 



o- = 2B4 



aE^ 



(6.25) 



This system of linear ordinary differential equations has 
the general solution 



and — ~A . With only the Hamiltonian (unitary) 
contributions, the matrix E_ appearing in the general so- 
lution is anti-Hermitian. 

In contrast, the contributions from the damping terms 



ait) 



e= (a(0)-a°)+aO 



E+t 

e= 



e= (a(0) 



(6.26) 



where gP satisfies E_qP — —A^^"^ , and the skew-symmetric 

matrix g° satisfies Rg° + g? = -2^- Note that if ^ 
is Hermitian and negative definite, then the dynamics 
will consist of some initial transients with a decay to the 
steady state: a(oo) = a°, g_{oo) = 

The first- and second-order physical moments also have 
a simple analytic form: 

{a){t) = e= ((a)(0) -a°)+a° 



e= ((:aa^ :) (0) 



-F{0))e=' + Fit). 



.27) 



A. 



=K 



r,(2)Q(2). 

^K ^K 



_r,(i)r,(2)* 

^K ^K 



^K ^K '-'ii' ^ T, 



=K 



^K ^K 
^K '^K 



'K "-^K 



..32) 



and C = —A — B . With only these damping con- 

==K =K o 
tributions, the matrix E is Hermitian. 



D. Bogoliubov dynamics 

As an example of how the dynamics of a linear problem 
can be solved exactly with these methods, consider the 



19 



following quadratic Hamiltonian 



(6.33) 



where x is a complex symmetric matrix. In the single- 
mode case, this Hamiltonian describes two-photon down- 
conversion from an undepleted (classical) pump|38i|. The 
full multi-mode model describes quasi-particle excita- 
tion in a BEC within the Bogoliubov approximation's^. 
Alternatively, it may be used to describe the dissocia- 
tion of a large molecular condensate into its constituent 
atoms Recasting this system into the general master- 
equation form [Eq. (|6.14() ]. we find that the constant and 
linear terms vanish, C — and 



where u) = lj — 17^/2. The block-diagonal form of these 
matrices allows us to write the solution to the phase-space 
equations as 



a{t) 


= e-''^*a(0). 


a+{t) 


= a+(0)e''^'*, 


nit) 




m{t) 




m+{t) 


= e''^^^*m+(0)e*"^* 



(6.39) 



If 7 is positive definite and commutes with oj, then the 
dynamics will be transient, and all these moments will 
decay to zero. 



E = -2B = 2A = 



X* 
X 



(6.34) 



The general solution [Eq. 



can then be written 



a{t) 



cosh*(|x|i) X 



■ sinh(|x|t) 



cosh* {\x\t) X 

sinh(|x|t) 
h*( 



cosh(|x|t) 



» sinh(|x|*) 



X cosh(|x|t) 

cosh*(|xK) X^^^^iJ^ 
X cosh(|x|t) 



a(0), 
g(0) 



(6.35) 



where the matrix cosh and sinh functions are as defined 
in Eq. (|5.26|l . If the system starts in the vacuum, for 
example, then the first-order moments will remain zero, 
whereas the second-order moments will grow as 



(:aSt:) = icosh*(2|xW-^^, 
(aa ) = -X rn ■ 



(6.36) 



E. Dynamics of a Bose gas in a lossy trap 

As a second example, we consider a trapped, noninter- 
acting Bose gas with loss modeled by an inhomogeneous 
coupling to a zero-temperature reservoir [37): 



d ^ 



LOijalaj,p 



-^7y (Saipa) ' 



(6.37) 

where a; is an Hermitian matrix that describes the mode 
couplings and frequencies of the isolated system, and 7 is 
an Hermitian matrix that describes the inhomogeneous 
atom loss. Recasting this in the general form, we find 
that = Tr7, A = and 



= 2 



T 

7 





, C = -i 


UJ 








7_ 







— U3* 



(6.38) 



F. Parametric amplifier 

A single-mode example that includes features of the 
previous two systems is a parametric amplifier consisting 
of a single cavity mode parametrically pumped (at rate x) 
via down-conversion of a classical input field and subject 
to one-photon loss (at rate ^'ifjlj: 

d \ 1 

'OlP"' 2 [^"^"^ ~ ^*°'°'' ^] + 2^ {2apa) - a)ap- pa) a) . 

(6.40) 

This corresponds to phase-space equations with 



= 2 



X* 
X 



= 2 



7 

-X 



7 



giving the solutions, for real x, 

a(t) = 



1 
1 
■(6.41) 



-7t/2 



-7t/2 



cosh xi sinh x* 
sinh xi cosh xt 



cosh x^ sinh x^ 
sinh xt cosh xt 



a(0) 



cosh x^ sinh xi 
sinh xt cosh x^ 



7 - '^X 



7^-2x^ 
XI J 



XI 
2 2x2 



(6.42) 



which are valid for 7 ^ 2x- For 7 > 2x, the system 
reaches the steady state a = 0, a; = a;*^, i.e. (a)a}j — 

■nP = 2x^/(7^ - 4x2) and (aa)° = = X7/(7^ - 4x^)- 

While this result is well-known and can be obtained 
in other wavs|i41,. ..42j. it is important to understand the 
significance of the result in terms of phase-space distribu- 
tions. In all previous approaches to this problem using 
phase-space techniques, the dynamically changing vari- 
ances meant that all distributions would necessarily have 
a finite width and thus a finite sampling error. However, 
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the Gaussian phase-space representation is able to han- 
dle all the linear terms in the master equation simply 
by adjusting the variance of the basis set. This implies 
that there is no sampling error in a numerical simulation 
of this problem. Sampling error can only occur if there 
are nonlinear terms in the master equation. These is- 
sues relating to nonlinear evolution will be treated in a 
subsequent publication. 



VII. CONCLUSION 

The operator representations introduced here represent 
the largest class of bosonic representations that can be 
constructed using an operator basis that is Gaussian in 
the elementary annihilation and creation operators. In 
this sense, they give an appropriate generalization to the 
phase-space methods that started with the Wigner rep- 
resentation. There are a number of advantages inherent 
in this enlarged class. 

Since the basis set is now very adaptable, it allows a 
closer match between the physical density matrix and 
appropriately chosen members of the basis. This implies 
that it should generally be feasible to have a relatively 
much narrower distribution over the basis set for any 
given density matrix. Thus, there can be great practi- 
cal advantages in using this type of basis for computer 
simulations. Sampling errors typically scale as l/^/T 
for an ensemble of T trajectories, so reducing the sam- 
pling error gives potentially a quadratic improvement in 
the simulation time through reduction in the ensemble 
size. As many-body simulations are extremely computer- 
intensive, both in real and imaginary time, this could pro- 
vide substantial improvements. Given the currently pro- 
jected limitations on computer hardware performance, 
improvement through basis refinement may prove essen- 
tial in practical simulations. 

We have derived the identities which are essential for 
first-principles calculations of the time evolution of quan- 
tum systems, both dynamical (real time) and canonical 
(imaginary time). Any quadratic master equation has 
an exact solution than can be written down immediately 
from the general form that we have derived. Higher or- 
der problems with nonlinear time-evolution can be solved 
by use of stochastic sampling methods, since we have 
shown that all Hamiltonians up to quartic order can be 
transformed into a second-order Fokker-Planck equation, 
provided a suitable gauge is chosen that eliminates all 
boundary terms. Because the Gaussian basis is ana- 
lytic, methods previously used for the stochastic gauge 
positive-P representation arc therefore applicable for the 
development of a positive semi-definite diffusion and cor- 
responding stochastic equations^, E| here. The abil- 
ity to potentially transform all possible Hamiltonians of 
quartic order into stochastic equations did not exist in 
previous representations. 

However, we can point already to a clear advantage to 
the present method in terms of deterministic evolution. 



For example, the initial condition and complete time evo- 
lution of either a squeezed state (linear evolution in real 
time) or a thermal state (linear evolution in imaginary 
time), with a quadratic master equation, is totally de- 
terministic with the present method. By comparison, 
any previously used technique would result in stochas- 
tic equations or stochastic initial conditions, with a finite 
sampling error, in either case. While this is not an issue 
when treating problems with a known analytic solution, 
it means that in more demanding problems it is possible 
to develop simulation techniques in which the quadratic 
terms only give rise to deterministic rather than random 
contributions to the simulation, thus removing the corre- 
sponding sampling errors. 

Finally, we note that the generality of the Gaussian 
formalism opens up the possibility of extending these rep- 
resentations to fermionic systems. 
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Appendix A: BOSONIC IDENTITIES 

To obtain the operator identities required to treat the 
time evolution of a general Gaussian operator, we need a 
set of theorems and results about operator commutators. 
These can then be used to obtain the result of the action 
of any given quadratic operator on any Gaussian opera- 
tor, as described in the main text. We use the following 
bosonic identities which are known in the literature, but 
reproduced here for ease of reference: 

Commutation Theorem (I): 

Given an analytic function p(a) with a power series 
expansion valid everywhere, the following commu- 
tation rules hold: 



p(a), 
ar, p(a^)] 



dp{a) 

dai 
dp{sL^) 

da 



,t 



(Al) 



Proof: Using a Taylor series expansion of p(a) 
around the origin in o^, one can evaluate the com- 
mutator of each term in the power series. Hence, 



p(a), ai^ 



E 

n 

E 

n>0 

dp{a.) 
dai 



rj I 



Pn 

7 a,- 



(n-1) 



(A2) 
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The second result follows by taking the Hcrmitian 
conjugate. 

Ordering Theorem (II): 

Given any analytic normally ordered operator func- 
tion p(a^ , a) with a power series expansion, the fol- 
lowing ordering rules hold: 



p(St,S)(al)" = 
(ai)Xa^a) = p(at,a) 



d_ 

da. 



p(a^a) 



(A3) 



Here the left arrow of the differential operator in- 
dicates the direction of differentiation. 
We can write these two identities in a unified form 
by introducing an antinormal ordering bracket, de- 
noted {: p : a}, which places all operators in anti- 
normal order relative to the normal term :p:. With 
this notation, we can write a single ordering rule for 
all cases: 



{: p(a) : a^} 



d 
da}, 



(A4) 



Proof: Since Oi commutes with all other annihilation op- 
erators and commutes with creation operators. 
Theorem (I) also holds for any normally ordered 
operator p(a^,a), with a power series expansion; 
provided derivatives are interpreted as normally or- 
dered also. The first case above then follows di- 
rectly from Theorem (I): 



t _ 



_d_ 

dai 



p(at,a) 



(A5) 



The required result then follows by using the equa- 
tion above n times, recursively. The second result 
is the Hermitian conjugate of the first. The last re- 
sult, Eq (|A4I) is simply a unified form that recreates 
the previous two equations. This can be applied re- 
cursively, since the RHS of this equation is always 
normally ordered by construction. 

Corollary: The antinormal combination of a Gaussian 
operator Ag (a) and any single creation or annihila- 
tion operator is given by a direct application of the 
ordering theorem, Eq HA3|I : 



|: Ag(a) : flpj =: [a^ 



A, 



(A6) 



It should be noticed here that the above expression 
assumes the covariance has the usual symmetry: 
then every operator occurs twice in the Gaussian 
quadratic term, which cancels the factor of two in 
the exponent. 



In the main text, these results are used directly to ob- 
tain all the required operator identities on the Gaussian 
operators. 



Appendix B: GAUSSIAN INTEGRALS 

In deriving the normalization, moments and operator 
identities of the Gaussian representations, we have had to 
calculate nonstandard integrals of complex, multidimen- 
sional Gaussian functions. The basic Gaussian integral 
that must be evaluated is of the form 



I = d 



-72 A/, 



(Bl) 



where, as in Eq. H3.2I) . the covariance ct is a 2M x 2M 
non-Hermitian matrix, and dz and Sz'^ are complex vec- 
tors of length 2M. There are two major differences be- 
tween this expression and the better known form of the 
Gaussian integral. First, the vectors dz = z — a and 
Sz^ = z* — contain offsets which are not complex- 
conjugate: a* . Second, the vector z does not 
consist of 2M independent complex numbers. Rather, 
it contains M independent complex numbers z and their 
conjugates z*. 

To evaluate such an integral, we first write it explicitly 
in terms of real variables as 



j2M 



ze 



where x_ — L_z_ — (Rez,Imz), = La = 

((a + a+)/2, (a - Q:+)/2j), and x = ' 
transformation matrix 



(B3) 



Note that the offset vector Xq will be complex, unless 
OL* = OL^ . We may remove it by changing variables u = 
x — XjQ and using contour integration methods to convert 
the integral back into an integral on a real manifold. 

With the offset removed, the square of the integral can 
be written in the form of a standard multidimensional 
Gaussian: 




ye 



d'^'ue 



-x^ x/2 ~y'^T_~^y/2 



- u/2 



(B4) 



where u = x + iy- Assuming that the matrix t ^ can be 
diagonalized: A = Hz"^ U} , we can factor the integral 
into a product of 2M integrals over the complex plane: 



2M 

n 



d w^e 



iLWf^/2 



(B5) 
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where w = U_u. These integrals can be evaluated by a 
transformation to radial coordinates, giving 



all eigenvalues of t have positive real part. Finally, noting 



that T = \L 



= 2-2*f 1 0- 1, we find that 



ir, \2M I I 

(27r) T 



(B7) 



(B6) 



which holds provided that all the Re Xuu > 0, i. e. that 



with the condition that the eigenvalues of |a;| have posi- 
tive real part. 
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